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Abstract 

The linear finite element approximation of a general linear diffusion problem with arbitrary 
anisotropic meshes is considered. The conditioning of the resultant stiffness matrix and the Jacobi 
preconditioned stiffness matrix is investigated using a density function approach proposed by Fried 
in 1973. It is shown that the approach can be made mathematically rigorous for general domains and 
used to develop bounds on the smallest eigenvalue and the condition number that are sharper than 
existing estimates in one and two dimensions and comparable in three and higher dimensions. The 
new results reveal that the mesh concentration near the boundary has less influence on the condition 
number than the mesh concentration in the interior of the domain. This is especially true for the 
Jacobi preconditioned system where the former has little or almost no influence on the condition 
number. Numerical examples are presented. 

Keywords: conditioning, flnite element, anisotropic diffusion, anisotropic mesh, stiffness matrix, 
extreme eigenvalue, Jacobi preconditioning, diagonal scaling. 
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1 Introduction 



Mesh adaptation is a common tool for use in the numerical solution of partial differential equations 
(PDEs) to enhance computational efficiency. It often results in nonuniform meshes whose elements vary 
significantly in size and shape from place to place on the physical domain. Nonuniform meshes could lead 
to ill-conditioned linear systems and their solution may deteriorate the efficiency of the entire computation. 
It is thus important in practice as well as in theory to understand how mesh nonuniformity affects the 
conditioning of linear systems resulting from discretization of PDEs on nonuniform meshes. The issue has 
been studied by a number of researchers mostly for the linear finite element approximation of the Laplace 
operator or a general diffusion operator; e.g., see |BS89l |DWZ09l |Pri73l |KHX12[ |She02l |Wat87| . More 



specifically, sharp estimates are available for the largest eigenvalue of the stiffness matrix DWZ09 Pri73 



KHX12 She02 . Indeed, the largest eigenvalue is bounded by a multiple (with a constant depending on 
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mesh connectivity) of the maximum of the largest eigenvalue of the local stiffness matrices. Moreover, 
the largest diagonal entry of the stiffness matrix is a good estimate for the largest eigenvalue KHX12| : 
it is tight within a factor of ci + 1 for any mesh, where d is the dimension of the physical domain. 

On the other hand, the estimation of the smallest eigenvalue is more challenging. Currently there are 
two approaches for this purpose. The first one utilizes Sobolev's inequality and is first used by Bank 
and Scott |BS89 for the Laplace operator with isotropic meshes. They develop a lower bound on the 
smallest eigenvalue of an appropriately scaled stiffness matrix and show that the condition number 
of the scaled stiffness matrix is comparable to that with a uniform mesh. The approach is recently 
employed by Kamenski, Huang, and Xu KHX12 to study a general diffusion operator discretized on an 
arbitrary anisotropic mesh. They establish bounds for the condition numbers of the stiffness matrix and 
the Jacobi preconditioned stiffness matrix and show that those bounds depend on three factors, one 
representing the condition number of the linear finite element equations for the Laplace operator on a 
uniform mesh and the other factors arising from the nonuniformity of the mesh viewed in the metric 
specified by the inverse of the diffusion matrix D (mesh -nonuniformity) and the mesh nonuniformity 
in volume measured in the Euclidean metric {volume-nonuniform,ity). It is also shown in [KHX12j 
that the Jacobi preconditioning — an optimal diagonal scaling for a symmetric positive definite sparse 
matrix — eliminates the effect of the mesh volume-nonuniformity and reduces the effect of the mesh 
D~ ^-nonuniformity. 

In the second approach (referred to as the density function approach hereafter), a lower bound on the 
smallest eigenvalue of the stiffness matrix is obtained through a lower bound of the smallest eigenvalue of 
a continuous generalized eigenvalue problem involving an auxiliary density function. It is first employed 
by Pried Pri73 for the Laplace operator. Pried shows how to obtain the lower bound for the continuous 
problem for spherical domains and claims via a physical argument that the bound is also valid for general 
domains. The obtained bound for the smallest eigenvalue of the stiffness matrix is valid for any mesh 
in any dimension but it is less sharp than that obtained in |KHX12 in three and higher dimensions; 
Fried's bound is proportional to power (d — 2) /d of the volume of the smallest mesh element whereas 
that in [KHX12 is only proportional to an average of the volumes of all elements (cf. (31) and (35)). 

In this paper we use the density function approach to develop bounds on the smallest eigenvalue 
and the condition number of both the stiffness matrix and the Jacobi preconditioned stiffness matrix 
for the linear finite element approximation of a general diffusion problem with arbitrary nonuniform 
meshes. The new bound on the smallest eigenvalue of the stiffness matrix is much sharper than that in 
the original work by Pried [Pri73j. Moreover, the new results are sharper than those obtained in KHX12 



in one and two dimensions and comparable in three and higher dimensions. In particular, they involve a 
factor describing the distance of a given element to the boundary and reveal that the mesh concentration 
near the boundary has less influence on the condition number than the mesh concentration in the interior 
of the domain. This is especially true for the Jacobi preconditioned system where the former has little 
or almost no influence on the condition number. The study shows that the density function approach 
can be made mathematically rigorous for diffusion problems with a general diffusion matrix and general 
domains for which an analytical expression of Green's function is generally unavailable. 

The outline of the paper is as follows. The linear finite element approximation of the general diffusion 
problem is described in Sect. [2j The lower bound for the smallest eigenvalue of the stiffness matrix 
is developed in Sect. [3j followed by the derivation of the bounds on the condition number in Sect. |4| 
Section |4] also contains comparison of the obtained bounds to those in |Pri73) and ^KIIX12j. Numerical 
examples are presented in Sect. |5j The conclusions are drawn in Sect. [6} 
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2 Finite element approximation 



We consider the boundary value problem (BVP) of a prototype anisotropic diffusion equation in the form 

r-V-(DVu) = /, inn, 
= 0, on dn, 

where is a simply connected polygonal or polyhedral domain in M'^ {d > 1) and D = D(a;) is the 
diffusion matrix. We assume that D is symmetric and positive definite and there exist two positive 
constants dmin and dmax such that 

dminl < B{x) < d^^^I, yx G 17, (2) 

where I is the d x d identity matrix and the less-than-or-equal sign means that the difference between the 
right-hand side and left-hand side terms is positive semidefinite. The weak formulation of the BVP ([T| 
is to find u £ V = -frd(O) such that 

A{u,v) = J^{v), yveV, (3) 

where 

A{u,v) = / \7v-B>\7udx and J^{v) = / fvdx. 
Jn Jn 

Let an affine family of simplicial meshes {Th} for be given. For any element X G 7^, let Fx : K ^ K 
be the affine mapping from the reference element K to the mesh element K and F'j^ the Jacobian matrix 
of Fk- For notational simplicity, we assume that \K\ = 1. Note that F'j^ is constant on K. 

Denote the linear finite element space associated with mesh Th by C V . Then, a linear finite 
element solution Uh G of BVP ([T| is defined by 

A{uh,Vh)=Hvh), yvh&V''. (4) 

Let N and be the numbers of the elements and the interior vertices of Th, respectively. If we order 
the vertices in such a way that the first vertices are interior vertices and denote by the standard 
linear basis function associated with the j**^ vertex, then we can express and Uh as 

= span{ </>!,... ,(/)Ar^ J and Uh = ^Uj(j)j. (5) 

i=i 

In the following we will use the function form Uh = Uj(j)j and the vector form u = [ui, . . . , uj^^^ 
synonymously. Using this, we can write Q in a matrix form as 

Au = f, 

where the stiffness matrix A and the right-hand side term / are given by 

Aij = A{(l}j,(t)i) = J2 1^1 ^'^i • ^Ky<Pj, i, j = 1, . . . , Nyi, (6) 

Ken 



fi= J2 I f(t>idx, i = l,...,N^^, 

KeTh 



IK 

and Ox is the average of D(a;) over K, i.e.. 



= -i- / Wx)dx. 



\K\ 



K 



We are interested in the condition numbers of matrices A and S ^AS ^, where S = diag {^/Ajj) is the 
Jacobi preconditioner. 
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3 Smallest eigenvalue of the stiffness matrix 

In this section we develop lower bounds for Amin(^) and Xram{S~^ AS^^) using the density function 



approach proposed by Fried Fri73 , where S = diag {^yAjj) is the Jacobi preconditioner. The approach 

(7) 



utilizes an auxiliary eigenvalue problem 

j —Au = Xpu, in J7, 



[u = 0, on d^, 

where p is a density function (to be chosen) satisfying 

< 

Pmin ^ P ^ Pmax < and J p{x) dx — 1. (8) 

The corresponding Galerkin formulation reads as 

( Vn, Vv) = \{pu, v), yv £ Hq (n) . 

Let Xp be the smallest eigenvalue of this eigenvalue problem and Amin(-Bp) the smallest eigenvalue of the 
Galerkin mass matrix Bp associated with the density p. Then, for any vector u, 

u^Au = {BVuh,Vuh) 

> dmin(Vu/i, Vn/i) 

> dminAp {pUh, Uh) 

^ '^min'^pAjnin(i3p) ||w||2 , 

which leads to 

-^min(^) ^ '^min-^pAjjjiji(5p) . (9) 

Thus, the key to the approach is to estimate Ap and choose the density function so that the lower bound 
([9]) is as large as possible. 

Hereafter, we assume that p is piecewise constant, i.e., p\k = Pk = constant, € Th- We denote 

l^minlp = min \ujj\ = ^ \K\^ = pK \K\ . 

It is known |KHX12[ Sect. 3] that 

{d+l){d + 2 

Inserting this into (l9|) yields 



XmUBp) > , (10) 



3.1 Estimation of the smallest eigenvalue of the continuous eigenvalue problem 

We now derive a lower bound on Ap by means of Green's function. The derivation consists of three 



lemmas, with the bound being given in Lemma 3.3 



First, we need the property of the strict positivity (or negativity) of eigenfunctions associated with 
the smallest eigenvalue of ([7|. 
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Lemma 3.1. For any density distribution p satisfying ([s]), the smallest eigenvalue Xp of the eigenvalue 
problem ([T]) is simple and positive. Any corresponding eigenfunction is either strictly positive or strictly 
negative in 



Proof. The positiveness of Xp follows from 



(Vu, Vu) Cp ||u|L 

Art = mm — ^ — > mm 



Cp 



u^O {pu, u) M^O 



>0, 



where Cp is the constant associated with the Poincare's i nequali ty. The rest of the proof is the same as 
for the standard eigenvalue problem with p = 1 (e.g., see 

replaced by 



GTOl 



Theorem 8.38]), except that ^^^'J"-* is 



□ 



Lemma 3.2. Let G{x,^) be the Green's function for —A subject to a homogeneous Dirichlet boundary 
condition on 0, and d{x) the distance from x to the boundary dO,, i.e., d{x) = min |a; — y|. Then, 



< C{x,^) < C X <^ ln(l + 



for d = 1, 
for d = 2, 



(12) 



□ 



|^-^P"Vin{l,^M|)}, ford>3. 

Proof For example, see |ADN59t [DS04{ |GW82[ |IR05). 

Lemma 3.3. Assume that p is piecewise constant and p £ (^1, ^^^2). Let dx = maXxi^K d{x). Then the 
smallest eigenvalue Xp of the problem ^ is bounded from below by 



E\K\,dK 



K 



Xp>C X { 



l + j:\K\An^l + dKP ri 



K 



d-2 



d 

d+2p 



li^r^ d^~^' 



p-1 
P 



for d = 1, 
for d = 2, 

for d>3. 



(13) 



Proof. For a positive eigenfunction Up associated with Xp we have 

Upix) = Xp [ G{x,^)p{0up{0dC 
Jn 

Consider two cases: 1 < d < 2 and d > 3. 



(14) 



For 1 < d <2, applying the Cauchy-Schwarz inequality to (14), 
Up{x)<Xp( [ p{^)G{x,^fd^ 



Taking square, multiplying with p{x), and integrating both sides over 0, gives 



p{x)up{x) dx < X^ / p{x) 



p{OCix,0^d^ 



dx. I p[i)up(ifdt 

Jn 
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This yields 



Ap > / p{x) 



p{i)G{x,ifdi 



dx 



(15) 



In ID, using this and Lemma |3.2 
>^P>c([ p{x 



p{i)d{x)d{i) di 



dx 



-1 



Ci / p{x)d{x)dx] >c[Y^\K\^dK 



K 



This gives (13) for the ID case. 
In 2D, from Lemma 3.2 we have 



p{i)G{x^ifdi <C( p(01n2 fl + ^^^^"l di<c( KOln^ fl + ^^'l d^, 

\ \x - ^\ j Jn \ \x - ^\ J 

where hn denotes the diameter of Let B^{x) be the disk with center x and radius e. Then, using 
property ^ we have 

/ p{i)G{x, d^<C I p{i) In^ (l + d^ + C f 

Jn Jn\BJx) V e / Je 



in\Be{x) 
< Cln^ ( 1 + 



< Cln^ 1 + 



d{x)hn 

£2 

d{x)hQ 



+ Cpmax£ 



d{x)hn 



Pmax In^ 1 + 



d{x)hn 



dr 



1 + InM 1 + 



d{x)h^ 



Taking e = pmix, we get 



p{i)G{x,^fdi<C l + ln\l + d{x)hnpm..) 



Inserting this into (15), recalling from (|8]), and noticing that d{x) < dx for all x € K, we have 

Xp>Gll + J2\K\pln^l + dKhnPm..)] 



K 



This gives ( 13 ) for the 2D case upon absorbing in the generic constant G. 

For d > 3, we need a slightly different version of the estimate (15). In this case. Lemma 3.2 implies 
that G{x, •) G LP(J7) for any x e 0, and p G ( 1, j^j- Using Holder's inequality. 



Up{x)<\J / G{x,^rd^ / piO'u.iO'd^ 



where q satisfies p + ^ = 1- By multiplying with p{x), taking power q, and integrating both sides over 
O, we get 



/ p{xyu,{xrdx<xi f p{xr( f G{x,^rdA''dx f pi^ru^i^rd^ 

Jn Jn \Jn / Jn 



and therefore 



Ap> / pix^i / G{x,^rd^) dx 



(16) 
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With Lemma 3.2 we have, using a similar strategy as in 2D, 

/ G{x,^rd^<C I |^-^|(2-'^)^'min|l, "^^^^'^^ 
Jn Jq [ 



< C 



x-if-'^^p di + c I \x-i\ 

Jn\B^{x) 

<C r r"!- V(2""')P dr + Cd{xfhP e~'^P 
Jo 

'^-^'^'^^P + Cdixfh^e-'^P. 



(2-d)p 



< C 



£ 



1 



di 



d-{d-2)p 

Choosing e such that the terms on the right-hand side are equal, we get 



e=({d-{d- 2)p) d{xYhl 



1 

d+2p 



and 



G{x,iYd^ < C{d-{d-2)p)~^d{xY^^^, 



where has been absorbed into the generic constant C. Inserting this into (16) leads to 



/ r dq q(d-(d~2)p) 

\p>C{j p{xY{d- {d-2)py^d{x) dx 
>C{d-{d-2)p)i+^v[ p{xYd{x) d+^p dx 



a ( q(d-(d-2)p) \ q 

>C{d-{d-2)p)^.\^pl\K\dj,'^'^ j • 



Since q = p/{p — 1) and p is element- wise constant, this gives (13) for d > 3. 



□ 



3.2 Smallest eigenvalue of the stiffness matrix 

Having established a lower bound on Ap, we can now proceed with the estimation of Amin(^) and 
Amin('S'~^^5'~^). Combining (11) with Lemma 3.3 we have 



Y.\K\,dK 



K 



max ) 



K 



d 

d+2p 



p d-(d-2)p^ 



p-1 



E\K\r' \KY^ dl7' ^+''' 



for d = 1 , 
ford = 2, (17) 

for d>3. 



The density function p is arbitrary so far. The optimal p is such that the right-hand side term of (17) 
attains the maximum value. It is difficult, if not impossible, to find the optimal p in general. We follow 
Pried [Fri73 to choose p such that 



\K\^ = const., \/K G Th- 



(18) 
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This gives 



and 



K 



N\K\ 



\K\ 



p N 



A I I ^ \ IT \ — P™'" \ 

ana I'^minlp ^ Prain |^min|p — ^ ^ ^5 



(19) 



(20) 



where Pmin is the minimum number of elements in a mesh patch. Inserting this into (17), using 
\K\ = \n\ /N (the average element volume) and Pmax = l/(-^ |-^min|)) we arrive at the following estimate. 

Lemma 3.4. The smallest eigenvalue of the stiffness matrix A for the linear finite element approximation 
of BVP ([T]) is bounded from below by 

-1 



Amin(-4) > ^ X < 



K 



N2^y\K\) % 



d+2p 



1 p d-(d-2)p ' 



P-1 



/or ci = 1, 
/or d = 2, (21) 

/or d > 3. 



Remark 3.1 (Optimality of the density function). For d 7^ 2, the choice of p in (19) is at least 
mong piecewisf 

Amin(^) > Cd 

n 



quasi-optimal among piecewise constant density functions. For the ID case, from (17) we have 

Pmin |-^min| 



T.\K\,dK 



> Cd„ 



K 



T.\K\,dK 



Cdr, 



Pn 



K 



K 



Maximizing the bound on the right-hand side is equivalent to minimizing the sum V^^- ,}^^'', dx- Since 

l-n-minL 



> 1 and dx > for all K, the optimal choice is 



analysis for the d >3 case is similar. 



1, which is equivalent to (18). The 



For the scaled stiffness matrix we have the following result. 

Lemma 3.5. The smallest eigenvalue of the Jacobi preconditioned stiffness matrix S'^AS~^ is bounded 
from below by 



N 



hE\K\/3KdK 



K 



Nd 



K 



hi:\K\PK] [jj1:\K\Pk l+ln2(l + d^7^ 



K 

d / p p d-(d-2)p ' 

\Ard(p-i) K 



p-i 



where Pk and ^k, defined as 



Pk 
Ih 



1 



maxi^ f^K 
Y.\K\f3K 

K 



for d = 1, 
ford = 2, (22) 
for d >3, 

(23) 
(24) 



reflect the non-uniformity of the mesh viewed in the metric specified by D ^ (see also Remark 4-3). 
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Proof. As for the non-scaled case we have (cf. (19]) and (10)) 



\miniS-^BpS-^)>Cmm'-^ 

and therefore 



Xrnin{S-^AS-^) > Cd^inApinin^. (25) 

It is known |KHX12| that 

Ajj <C J2 1^1 • \\iF'Kr^^K{F'K)-^^ = Cd^in E 1^1 (^K, J = 1, . . . ,iV«, (26) 

where C is a constant depending only on K and the hnear basis functions on K. Thus, 

. \^j\p ^ C . Ekgu^j \K\pk 
mm > mm — . 

j Ajj dmin i L^Kduij \K\ PK 

Once again, we choose p to get rid of the minimum sigiij^ viz., 

mm — = const. 

This leads to 

p{x)\k 



E \K\I3k 

KdTh 



and 

. Y.K>^u,,\K\ pK 1 
mm 



KdTh 



The estimate (22) follows from this, inequality (25) and Lemma 3.3 □ 



4 Condition numbers of stiffness matrix and Jacobi preconditioned 
stiffness matrix 



We first quote an estimate on the maximum eigenvalues of A and 5 ^AS ^ from KHX12 . 

Lemma 4.1 ( |KHX12] ). The largest eigenvalues of the stiffness matrix A and preconditioned stiff- 
ness matrix S~^AS~^ with the Jacobi preconditioner S = diag (y^Ajj) for the linear finite element 
approximation of B VP ([T| are hounded by 

max Ajj < Aniax(^) < (d + 1) max Ajj < dmin v d{d + 3) max \K\ Pk, (27) 

l<Amax(5-i^5-i)<d + l. (28) 



^Although it is not quite clear if this choice is optimal, it is the best choice we could find so far. 
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Combining Lemma 4.1 with Lemmas 3.4 and 3.5 we obtain our main theorem. 



Theorem 4.1. Denote the average element volume by \K\ 



N 



and let p £ {1 



' d~2 



I. Then the condition 



numbers of A and S ^AS ^ for the linear finite element approximation of EVP ([T]) are bounded by 



( 1 



k{A) < CNd ( N~ max 



and 



K 



X < 



d I M 



p d-(d-2)p ' 
d+2p 



£-1 

P 



for d=l, 
ford = 2, 

for d>3, 
(29) 



772 I -^1 PxdK, 



K 



wE\K\pKj i^j^E\K\(3K[l + ln^l + dK7h 

/ p p d— (d — 2)p \ p 

1 ' ^ EI^I/?F'^^F'' "^'^ 



for d = 1, 
for d = 2, 

for d > 3, 



(30) 



where Pk o^nd defined in (23) and (24|. 



Remark 4.1. From the theorem we can see that both bounds for n{A) and k[S~'^ AS~^) contain the 
maximum distance dx from element K to the boundary of the domain, dx becomes smaller when K 
is closer to d^. As a consequence, elements close to the boundary have less influence on the bounds 
than those away from the boundary. This information is useful for problems having boundary layers for 
which adaptive meshes are typically dense near the boundary. 



Remark 4.2 (Comparison to KHX12| ). The bounds in (29) and (30) are sharper than (in the cases of 
d < 2) or comparable to (in the case of d > 3) those obtained in KHX12| [Theorem 5.2] using Sobolev's 
inequality, 



k{A) < cm X ( iV"^ max ^ \K\I3k)^' 



1, 

1 + ln 



1^1 



N \\K\) 



d-2 
2 



for d = 1, 
for d = 2, 

for d > 3. 



(31) 



and 



k{S-^AS~^) < CNd X < 



,wYk\K\/3k), ford=l, 
wJ:K\K\PK){l + \ln-fh\), ford = 2, 

2 

d\d 

for d>3. 



N 



EK\K\f3k 



(32) 



For d < 2, bounds (31 ) and l\32n follow directly from Theorem 4.1 if we replace dx by its largest possible 



value hn (the diameter of i7). For d > 3, if we replace dx with hn and take p close to d/(d — 2), the new 



bounds are very close to the bound (31) and (32). In this sense, they are comparable. 
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Remark 4.3 (©"^-uniform meshes). Like (31) and (32), bounds (29) and (30) involve three factors: 
the basic factor N'^/'^ which corresponds to the condition number of the stiffness matrix for the Laplace 
operator on a uniform mesh, the factor involving (3k which reflects the mesh non-uniformity in the metric 
specified by D~^, and the factor involving the element volume which represents the mesh non-uniformity 
in the Euclidean metric. It is not difficult to show that 



f 1 



tr 



N 



<PK<d 



K\ det(DK 



N\-d 



det ( {F'K)-'^}K{F'K) 



I 



N 



|i^|det(DA-)"^ 



1 I ) 

2 / 



where Uh = J2 \^\ det(Di^) 2 . For any uniform mesh in the metric specified by D~ we have (e.g., see 

K 

Huang and Russell |HR11| ) 
1 



d 



tiUF^ 



i,)-'BK{F^)-^ 



and thus. 



From (29) and (30), we get 



det((F;,)-^D;,(F;,)-^^^ 



<PK<d 



and 



N 



\K\det{DK) 



k{A) < C 



2 

<7h 



hY.dK-, 



N 



K 



(iVl^maxI) X < 



N Y 



K 



1 + InM 1 + 



dK\K\ 



and 



d I M 



i:\K\dK, 

K 



N\i\\K\) % 



P a-(d-2)p \ p 
d+2p \ 



k{S-'AS-') < C 



N\-d 



X < 



El^l l + ln2(l + di^) 
K 



p d-(d-2)p 



p-1 



If we replace by /ij^ in ( 34 ) , we get 



k{S-^AS-^) < C 

which gives exactly the same bound obtained in KHX12 



2 

N\d 



for d = 1, 
for d = 2, 

for d > 3, 



for d = 1, 
for d = 2, 

for d > 3. 



for ©^""^-uniform meshes. 



(33) 



(34) 
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Remark 4.4. As mentioned before, the approach we used here to estimate Amin(^) and Amin(5'~^^5'~^) 
was first proposed by Pried Pri73 . However, there is significant difference between our development and 
Fried's. First, Fried obtains a bound on X'^i^(—A) for balls where the analytical formula of the Green's 
function is available and claims using a physical intuition that the bound is also valid for other domains. 
Our derivation for the bound on X'^^^^—A) is mathematically rigorous, as seen in Lemmas 3.1-3.3 
Second, the lower bound obtained in Fried |Fri73 for Anim(^) can be expressed in the current notation 
as 

1, for d = 1, 

k\ 



Amin(^) ^ 



N 



1 +ln 



|Xmm| 
2 

d ^ 



1 



for d 



for d>3. 



(35) 



It is easy to see that our estimate (21) is sharper than (35). Third, Fried estimates the maximum 



eigenvalue of the stiffness matrix by the maximum of the maximum eigenvalues of the local stiffness 
matrices, with the latter being computed directly for the Laplace operator. On the other hand, our 
estimate on the maximum eigenvalue of the stiffness matrix in Lemma [4. 1| is not only more general, valid 
for all second order diffusion operators, but also very accurate. Finally, we would like to mention that 



Fried |Fri73 does not study preconditioning for the stiffness matrix. 



5 Numerical experiments 



The dependence of the conditioning on the diffusion matrix was extensively discussed in [KHX12 and 
the purpose of this paper is to investigate the dependence of the conditioning on the mesh density 
throughout the domain. Thus, for the simplicity, we set D = / (Laplace operator) in all of our numerical 
examples. 

Bounds on the smallest eigenvalue contain a constant C independent of the mesh but dependent on the 
dimension, domain, reference element, and basis functions on the reference element. In our computation 
we obtain its value empirically by calibrating the bounds for Amm('S'~^^5'~^) through comparing the 
exact and estimated values for a series of uniform meshes. 



5.1 ID 



In ID, the new bounds (29) and (pO[) (with D = /) become 



k{A) < C ^ dx ■ max ^ 



K 



k{S-^AS-^) < CJ2 



K 



\K\ 



On the other hand, the bounds (31) and (32) from |KHX12| reduce to 

k{A) < CN ■ max ^ ^ 



i{S-^AS-^) < CJ2 



K 



1 



(36) 
(37) 

(38) 
(39) 



We consider two set of meshes formed by Chebyshev nodes and nodes simulating extreme boundary 
layer refinement. 
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5.1.1 Chebyshev nodes 



The first ID example is a non- uniform mesh given by the end points xq = and xn = 1 and the 
Chebyshev nodes, 



1 



= -(1 -cos(^j)), i = 1, . . . , - 1, 



(40) 



where 



Co = 0, Cat = 1 and Cj 



vr(2j - 1) 
2(iV - 1) 



j = l,...,iV-l. 



For this mesh, we have 



\Kj\ = Xj - Xj-i ~ ^sin(^j), dxj = min Q(l - cos(Cj)), ^(1 + cos(Cj_i))^ . 



Using these, we get 



and 



max 

j 



1 

\K\ 



T^TTT ~ max 



1 2 "^Z' 1 

K j=l 



(1 -cos(O) ~ 1. 



Thus, from (36) and (38) we see that the new bound and the bound from [KHX12 both lead to 

k{A) < CN^. 

For the scaled case, we have 

■f 1 — cos(C) 



E 

K 



\K\ 



fsin(e,) 



sin(^) 



and 



1 



1 



N/2 



N 



Thus, we have 



for the new bound ( 37 ) and 



Ci sin(^) 
k{S-^AS-^) < CN^ 
k{S-^AS-^) < CN^ In N 



~ N^ln\tan{^/2)\ ||^ ~ iV^lnA^. 



(41) 



from the bound (39) from KHX12 . Notice that (41) has the same order as k,{S~^AS~^) for a uniform 
mesh as N increases. 

The numerical comparison of the estimated and exact values for k(A) and k(S^^AS^^) are presented 
in Fig. [T| As expected, for the non-scaled case both estimates are comparable and very tight (Fig. la). 
On the other hand, the new estimate (29) is more accurate than the one from KHX12 if diagonal scaling 
is applied (Fig. lb). For this example, the new bound on k(S^^ AS~^) seems to be asymptotically exact. 
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— 1 — exact 










A new 










-■- |KHX12| 








uniform mesh 












^^^^ 















10^ 10^ 
(a) k{A) vs. number of elements N 



10^ 10^ 
(b) k{S^^ AS^^) vs. number of elements N 



Figure 1: Conditioning for the Chebyshev nodes (40) 
5.1.2 Boundary layer refinement 

The second ID example simulates boundary layer mesh refinement towards the boundary point x = 
with the internal nodes 

j = l,...,N-l. (42) 



2V2 



N 



For this mesh, we have 



2^ 



-1 



27V 



dK, =mm U^,l 



2^ 



-1 



2N 



Then, bounds (38) and (39) for KHX12 become 



k{A) < CN2^ and niS'^AS-^) < C2 



N 



whereas the new bounds ( |36| ) and (37) give 

k{A) < C2^ and k{S-^ AS'^) < CN. 

This shows that, for boundary layer refinement, the new bound for the scaled case is a significant 
improvement. Note that n^S""^ AS~^) = 0{N), which has a smaller order as — t- oo than the condition 
number for a uniform mesh (which is 0{N'^)). Thus, for ID problems with steep boundary layers, strong 
mesh concentration towards the boundary not only increases the accuracy of the solution but, at the 
same time, improves the conditioning. 



Numerical results are in perfect agreement with the analysis. Indeed, Fig. 2a shows that both new 
and old estimates are comparable for the non-scaled case, with the new one being slightly more accurate. 



After scaling (Fig. 2b), the situation is quite different: the new estimate is very close to the exact value 



whereas the estimate from |KHX12 exhibits a dramatic overestimation. 



5.2 2D 

In 2D we consider a mesh for the unit square [0, 1] x [0, 1] with 0{^fN) skew elements near the boundary. 



as shown in Fig. 3a 
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3 5 9 17 

(a) k{A) vs. number of elements N 




3 5 9 17 33 

(b) k{S~^ AS~^) vs. number of elements N 



Figure 2: Conditioning for the boundary layer refinement (42) 



First, we set the maximum aspect ratio at 125 : 1 and verify the dependence of the condition number 
on the number of mesh elements N (Fig. [4]). Then, we set = 20, 000 and change the maximum aspect 
ratio of the mesh elements to investigate the dependence of the conditioning on the mesh shape (Fig. |5]) . 

The observation is that the new and the [KHX12 estimates are comparable, with the new one being 
shghtly more accurate, especially when changing the maximal mesh aspect ratio. However, neither of 
the estimates is as accurate as in ID, meaning that our estimation of the smallest eigenvalue can be 
further improved. This observation is essentially valid in higher dimensions as well. Since the Green's 
function has a singularity for d >2 (and depends on the shape of the domain) , it is difficult to obtain 
estimates as sharp as in ID in general. Interestingly, in our example, the exact condition number of the 
scaled stiffness matrix appears to be independent of the aspect ratio of the boundary layer elements. 
This suggests that the optimal estimate in 2D, similar to the ID case, should be 



k{S-'AS-') <CY,\K\f^K log (l + ^) 



(43) 



5.3 3D 



Similarly to the 2D case, we use a mesh for the unit cube [0, 1]'^ with 0(A^^/^) skew elements near the 
boundary (Fig. 3b) and consider two different settings: fixed anisotropy (25 : 25 : 1) with increasing 
number of elements and fixed N = 29, 478 paired with the changing anisotropy of the mesh. Numerical 
results for p = 2.9 are presented in Figs. [6] and [7} 



First, we observe that the new estimate is comparable to the estimate from KHX12 for both non-scaled 



and scaled cases. Second, as in the 2D case, we observe that increasing the maximum aspect ratio of 
elements at the boundary has no impact on the exact condition number of the Jacobi preconditioned 
stiffness matrix (at least for the considered mesh type). This indicates that the obtained bounds for 
the scaled stiffness matrix are not optimal. The numerical results also suggest that the optimal bound 
should have a much stronger dependence on the distance dx of an element K from the boundary of the 
domain. 
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(a) 2D meshes used in Sect. 



5.2 



(b) 3D meshes used in Sect. 



5.3 



Figure 3: Test meshes in two and three dimensions 
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(a) k{A) vs. number of elements N 
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(b) k{S^^ AS^^) vs. number of elements N 



Figure 6: 3D meshes with a fixed aspect ratio of 25 : 25 : 1 and changing N 
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10° 10^ 10^ 10^ 10° 10^ 10^ 10^ 

(a) k{A) vs. number of elements N (b) k{S~^ AS~^) vs. number of elements N 

Figure 7: 3D meshes with a fixed N = 29,478 and changing aspect ratio 



6 Conclusions 

In the previous sections we have studied the linear finite element approximation of the boundary value 
problem ([T]) with general nonuniform meshes and developed bounds for the condition numbers of the 
stiffness matrix and Jacobi preconditioned stiffness matrix. The main result (Theorem 4.1) shows that 
the density function approach of Fried [Fri73j can be made mathematically rigorous for general domains 
and lead to estimates that provide more detail and are sharper than existing estimates in one and two 



dimensions and comparable in three and higher dimensions. Moreover, the bounds in Theorem 4.1 
involves a factor dx which describes the maximum distance from element K to the boundary of the 
domain and becomes smaller when K is closer to di}. They reveal that the mesh concentration near 
the boundary has less influence on the condition number than the mesh concentration in the interior of 
the domain. This is especially true for the Jacobi preconditioned system where the former has little 
or almost no effect on the condition number. The numerical results presented in Sect. [5] confirm the 
theoretical analysis although they also suggest that the new bounds could be further improved in two 
and higher dimensions. 
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